/**
 * @file FEMFunction.h
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2021-01-17
 * 
 * @copyright Copyright (c) 2021
 * 
 */

#ifndef _FEMFUNCTION_H_
#define _FEMFUNCTION_H_
#include "Element.h"
#include "Point.h"
#include "Mesh.h"
#include "Dofs.h"
#include "Eigen/Sparse"
#include "Eigen/Dense"

#define TEMPLATE template<Dimension DIM>

TEMPLATE
/** \ingroup The FEMFunction for FEMSpace
    * \brief The FEMFunction for FEMSpace.
    * 
    *      This class is designed to compute the 
    *       u_h = sum_{i = 0}^{n}u_i * phi_{i}.
    *       for one point in Dimensional DIM p,
    *      get the value u_h(p) or the gradient grad(u_h)(p) or other scalar and vector;  
    * \tparam DIM The Dimensional of this function.
 */
class FEMFunction
{
public:
    using VectorXd = Eigen::VectorXd;
    using Element = Element<DIM>;
    using Element_pointer = Element*;
    using const_Element_pointer = const Element *;
    using Element_ref = Element&;
    using const_Element_ref = const Element&;

public:
    FEMFunction(){};
    FEMFunction(Element_pointer element,VectorXd v):_element(element),_v(v){};
    FEMFunction(VectorXd v):_v(v),_element(NULL){};
    ~FEMFunction(){};
    void SetElement(Element_pointer element){_element = element;};

private:
    Element_pointer _element;
    VectorXd _v;
public:
    /**
    * @brief Compute the value of u_h at point p in interal element;
    * 
    * @param p Point<DIM> The point must be in the reference element;
    * @return double 
    */
    double value(Point<DIM> p);
    /**
     * @brief Compute the value of u_h at point p in the vector in interal element;
     * 
     * @param p_v The point in the vector must be in the reference element;
     * @return std::vector<double> 
     */
    std::vector<double> value(std::vector<Point<DIM> > p_v);

    std::vector<double> gradient(Point<DIM> p);

    std::vector<std::vector<double> > gradient(std::vector<Point<DIM>> p_v);

    double value(Point<DIM> p,Element_pointer element);

    std::vector<double> value(std::vector<Point<DIM> > p_v,Element_pointer element);

    std::vector<double> gradient(Point<DIM> p,Element_pointer element);

    std::vector<std::vector<double> > gradient(std::vector<Point<DIM>> p_v,Element_pointer element);

    double& operator()(int i)
    {
        if(i > _v.size())
            std::exit(-1);
        return _v[i];
    }

    double& operator[](int i)
    {
        if(i > _v.size())
            std::exit(-1);
        return _v[i];
    }
};
TEMPLATE
double FEMFunction<DIM>::value(Point<DIM> p)
{
    double ans = 0.0;
    if(DIM == 2)
    {
        double xi = p[0];
        double eta = p[1];
        for(int i = 1;i <= _element->n_Dofs();i++)
            ans += _v(_element->NdIdx(i)) * _element->phi(xi,eta,i);
    }
    return ans;
}
TEMPLATE
std::vector<double> FEMFunction<DIM>::value(std::vector<Point<DIM> > p_v)
{
    std::vector<double> ans_v;
    if(DIM == 2)
    {
        double xi = 0.0;
        double eta = 0.0;
        double ans = 0.0;
        for(int j = 0;j < p_v.size();j++)
        {
            xi = p_v[j][0];
            eta = p_v[j][1];
            for(int i = 1;i <= _element ->n_Dofs();i++)
                ans+=_v(_element->NdIdx(i)) * _element->phi(xi,eta,i);
            ans_v.push_back(ans);
        }
    }
    return ans_v;
}
TEMPLATE
std::vector<double> FEMFunction<DIM>::gradient(Point<DIM> p)
{
    std::vector<double> ans_v(DIM);
    if(DIM == 2)
    {
        double xi = p[0];
        double eta = p[1];
        double gradient_x = 0.0;
        double gradient_y = 0.0;
        for(int i = 1;i <= _element ->n_Dofs();i++)
        {
            gradient_x += _v(_element->NdIdx(i)) * _element -> phi_x(xi,eta,i);
            gradient_y += _v(_element->NdIdx(i)) * _element -> phi_y(xi,eta,i);
        }
        ans_v[0] = gradient_x;
        ans_v[1] = gradient_y;
    }
    return ans_v;
}
TEMPLATE
std::vector<std::vector<double> > FEMFunction<DIM>::gradient(std::vector<Point<DIM> > p_v)
{
    std::vector<std::vector<double> > ans_vv;
    ans_vv.resize(p_v.size());
    for(int k = 0;k < ans_vv.size();k++)
        ans_vv[k].resize(DIM);
    if(DIM == 2)
    {
        for(int j = 0;j < p_v.size();j++)
        {
            for(int i = 1;i <= _element -> n_Dofs();i++)
            {
                double xi = p_v[j][0];
                double eta = p_v[j][1];
                ans_vv[j][0] += _v(_element->NdIdx(i)) * _element -> phi_x(xi,eta,i);
                ans_vv[j][1] += _v(_element->NdIdx(i)) * _element -> phi_y(xi,eta,i);
            }
        }
    }
    return ans_vv;
}
TEMPLATE
double FEMFunction<DIM>::value(Point<DIM> p,Element_pointer element)
{
    double ans = 0.0;
    if(DIM == 2)
    {
        double xi = p[0];
        double eta = p[1];
        for(int i = 1;i <= element->n_Dofs();i++)
            ans += _v(element->NdIdx(i)) * element->phi(xi,eta,i);
    }
    return ans;
}
TEMPLATE
std::vector<double> FEMFunction<DIM>::value(std::vector<Point<DIM> > p_v,Element_pointer element)
{
    std::vector<double> ans_v;
    if(DIM == 2)
    {
        double xi = 0.0;
        double eta = 0.0;
        double ans = 0.0;
        for(int j = 0;j < p_v.size();j++)
        {
            xi = p_v[j][0];
            eta = p_v[j][1];
            for(int i = 1;i <= element ->n_Dofs();i++)
                ans+=_v(element->NdIdx(i)) * element->phi(xi,eta,i);
            ans_v.push_back(ans);
        }
    }
    return ans_v;
}
TEMPLATE
std::vector<double> FEMFunction<DIM>::gradient(Point<DIM> p,Element_pointer element)
{
    std::vector<double> ans_v(DIM);
    if(DIM == 2)
    {
        double xi = p[0];
        double eta = p[1];
        double gradient_x = 0.0;
        double gradient_y = 0.0;
        for(int i = 1;i <= element ->n_Dofs();i++)
        {
            gradient_x += _v(element->NdIdx(i)) * element -> phi_x(xi,eta,i);
            gradient_y += _v(element->NdIdx(i)) * element -> phi_y(xi,eta,i);
        }
        ans_v[0] = gradient_x;
        ans_v[1] = gradient_y;
    }
    return ans_v;
}
TEMPLATE
std::vector<std::vector<double> > FEMFunction<DIM>::gradient(std::vector<Point<DIM>> p_v,Element_pointer element)
{
    std::vector<std::vector<double> > ans_vv;
    ans_vv.resize(p_v.size());
    for(int k = 0;k < ans_vv.size();k++)
        ans_vv[k].resize(DIM);
    if(DIM == 2)
    {
        for(int j = 0;j < p_v.size();j++)
        {
            for(int i = 1;i <= element -> n_Dofs();i++)
            {
                double xi = p_v[j][0];
                double eta = p_v[j][1];
                ans_vv[j][0] += v(_element->NdIdx(i)) * element -> phi_x(xi,eta,i);
                ans_vv[j][1] += v(_element->NdIdx(i)) * element -> phi_y(xi,eta,i);
            }
        }
    }
    return ans_vv;
}

#else
#endif